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INTRODUCTION 


Solutions of constrained optimization problems minimize an objective 
function, E, subject to given constraints. In aerodynamic applications, the 
objective function and the constraint functions, fj , i = 1 , 2, . . ., depend 
on the flow field solution, g. The optimization scheme developed here is 
applicable to situations in which the flow governing equations are nonlinear 
equations that are solved iteratively. 

Conventional optimization methods (e.g., the steepest descent method and 

the conjugate gradient method) are iterative procedures that require the 

evaluation of the objective function many times before the converged optimum 

solution is determined. Since E and f are dependent on the flow solution, g, 

in addition to the vector of design parameters, P, the flow governing equation 

must be solved each time E and f are evaluated. Therefore, the application of 

1 1-5 

conventional optimization schemes to aerodynamic design problems leads to 
two-cycle (inner-outer) iterative procedures. The inner iterative cycle solves 
the analysis problem for g iteratively, while the outer cycle determines the 
the optimum P iteratively. An alternative to this costly procedure is the 
approach based on the idea of updating the flow variable iterative solutions 
and the design parameter iterative solutions simultaneously. 

In this manual we present a description of the optimization scheme. The 
formulation presented is that for optimizing the design of a propeller 
by maximizing its efficiency, subject to a given power constraint. The Euler 
equations are assumed to be the flow governing equations. An implicit 
approximate factorization scheme is used to compute the flow field about 
the propeller. 
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TECHNICAL BACKGROUND 


The propeller design problen is cast into an optimization formulation in 

* 

which the optimum design parameter vector, P , is to be determined such that 


subject to the constraint 


E C P* ; g) = rain E(P; g) 
P 


f(P; g) = 0 


( 1 ) 

( 2 ) 


with the flow variable vector g satisfying the flow governing equation 


D(g) = 0 


(3) 


subject to the boundary condition 

B(g; P) - 0 (4) 

Our objective is to maximize the propeller efficiency, q . The objective 
function is therefore defined by 

E - -n 

The propeller power requirements are constrained to a specified value through 
the constraint function 


4 


4 


4 


4 


( 


f - C - C 


P Po 

Equation (3) is the system of Euler equations governing the flow field, and 
Equation (4) is the propeller solid wall boundary condition. The vector of 
design parameters, P, defines the propeller geometrical configuration. 

The goal of the optimization scheme is to determine the values of the 
design parameters that minimize the objective function, E, subject to an 
equality constraint. A search must therefore be conducted in the design 
parameter space P for the optimum solution, P*. This optimization problem is 

A 

most conveniently solved in the rotated design parameter space P, with the 
P^ coordinate normal to the constraint surface and the P^ coordinates, where 
S, = 2, 3, . . . , L, parallel to the constraint surface. For fixed values of 
the components of f, let 


n+1 

§ 


- !fc(g n ; p), n = o. 


1 , 2 , 


(5) 
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be the iterative solution for the analysis problem, where j denotes the 
solution obtained by applying the iterative scheme for solving the Euler 
equations once using g n as an initial guess. An implicit approximate factor- 
ization scheme is used here to solve the Euler equations. It is described in 
Reference 6. As for the analysis solution, obtaining the optimization solu- 
tion requires the repeated application of Equation (5) to update the flow 
field. While P is held fixed in the former case, it is allowed to vary in the 

a 

latter. The scheme used to update P follows. 

A 

The vector of design parameters P is updated every aN iterations. 
Therefore, 


where 


*n+l ~n 
P * P 


5P n+1 


( 6 ) 


6P n+1 - 0, (n+l)/AN + 1, 2, 3, . . . 


In the iterative steps that satisfy the relation (n+l)/AN - 1, 2, 3, . 
incremental values for the design parameters are given by 


where 


«p" +i 


«p" +1 - rain 11, 


In /l 


[min 

cc if i, 5P : 

If n | 

max 


max \ 

ap" , a = 2, : 

1*0/ 


t n - fO"; s"> 


-n+l 1 r f n+1 . - . ( n+1 1N , *n+l-AN 

4P, +D+C 2 (T t -D]«P t 


n+1 _ I ft 

‘ ’ ’ llEj 6f; +1 ' 4 "l 


the 


(7) 

( 8 ) 


(9) 


AE 


n 

l 


E(| n + e 



Here, e is a small positive constant and a m 1, 2, . . . , L, are the set of 

orthogonal unit vectors along the axes of the rotated coordinate system 

P~, . . P™. The solution is a solution in which the component 

1 Z L ft 

of P is perturbed by e. 
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The incremental displacement in the design parameter space, introduced so 
that the constraint nay be satisfied, is taken in the direction normal to the 
constraint surface and is determined by the chord method in Equation (7)* The 
constant sets an upper limit on the magnitude of this incremental dis- 

placement. The incremental displacements given by Equation (8) are introduced 
along the coordinate axes, which are parallel to the constraint surface with 
the purpose of reducing the value of the objective function. The sign of the 

incremental correction 6P*) + \ where 6P n+ ^ is the component of the vector 

* n+ 1 ^ ^ a n 

6P , is chosen to be opposite to that of aE/aP^. The magnitude of the incre- 
ment is given by 


6P" +1 

“ c 

5 £n+l-AN 

l 


* l 


with an upper limit given by <5P , where c > 0. If the signs of 6P™ + '*' and 

*n+l-AN n3X ~n 1 ^n-A’I 

5P^ are in agreement, then the last two iterative solutions P^ and P^ 

a 

fall to one side of the point along the P £ direction at which E is a minimum. 
In this case, c is set equal to the constant , which is greater than 1. 
Increasing the magnitude of the step size in this manner accelerates the 
approach toward the point along the P^ direction at which E is a minimum. On 
the other hand, if the signs of <$P^ + ^ and 5P™ + ^ are not in agreement, then 

p” and p! 1 j fall on opposite sides of the point along the P direction at 

K K, 

which E is a minimum. In this case, c is set equal to the constant c^, which 

is less than 1. Decreasing the magnitude of the step size in this manner is 

necessary for convergence to the point along the P^ direction at which E is a 
minimum. 

The updated components of the design parameter vector P are used to 
calculate the new flow iterative solution, g n+ \ given by 


n+1 , n ^n+1* 

g = jtjCs ; g ) 


( 10 ) 


and the perturbed solutions g!* + \ l “ 1, 2, . . . , L, given by 

— X, 


n+1 _ , / n ~n+l $n+l\ 
g„ - i(8„ ; P + e i„ ) 


(ID 


While the optimization procedure is most suitably conducted in terns of 

1, 2, . . . , L, the flow solution is 

1, 2, • • • , L. 


the transformed parameters P , i 


computed in terras of the physical design parameters P , g, 
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In order to express the transformed design parameters In Equations (10) and 
(11) In terms of the original design parameters. It Is necessary to use the 
transformation equation, which relates these two sets of parameters. This 
equation Is 

p n+l = T n+1 -n+1 

where the orthogonal transformation matrix T n+ ^ Is given by 


„n+l _ r . n+1 , n+1 


,n+l. 


U lirj. urx . U i X «» 

1 ~2 * * ~L J 

The unit vector i p is normal to the constraint surface at P = P and is given 
by 


n+1 

~1 


7fO* n ; g n ) 

IvfCfj g n )l 


(12) 


where an estimate for G the a th component of vf, is given by 


K 


ftp + t I"; s") - f(t"; ? n ) 

a ~a ~ 


(13) 


The Gram-Schmidt orthogonalization process, which uses a set of L linearly 
independent vectors to construct a set of L orthonorraal vectors, is used to 

Yl+1 

construct the unit vectors i , a “ 2, 3, . . . , L, along the rotated axes 
n+1 ^ 

P , A a 2, 3, . . I. The following equation is used for this purpose: 


.n+1 


where 


I n+1 

l^T’ 


a-l 


a - 2 , 3 , . . . , L 


if 1 - 1? - H Cl? • dp 1 ) C 1 


k-1 


In the initial iterative step, the vectors i are given by i “ e , 

a a a 

a ■ 1, 2, . . . , L, where e , a “ 1 , 2, . . . , L, are the set of orthogonal 
unit vectors along the axes of the coordinate system p , P , . . ., P . 

1 Z 

While the flow variable vector g is updated each iterative step, the 
coordinate system in the design parameter space is rotated every aN iterations. 
The unit vectors i , like the vector of design parameters P, are updated only 
in the iterative steps that satisfy the relation (n+l)/AN = 1, 2, 3, . . . 
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The optimization scheme described above requires that L+l iterative 
problems be solved in parallel. In addition to the main solution, L perturbed 
solutions are computed in which each of the design parameters in the trans- 
formed space P^ , § 2 * . . . , is perturbed. The computational costs and the 
computer memory requirements are therefore proportional to L+l. A modifica- 
tion to this scheme, which requires that only L iterative solutions be 
obtained, is now introduced. In the modified procedure, the perturbation 
solution associated with the perturbed design parameter in the direction of 
the P. axis, normal to the constraint surface, is not computed. This solution 

was used in Equation (13) to compute G-. , which is required for the calculation 
n+1 1 

of the vector i^ 3 which determines the direction normal to the constraint 
surface in Equation (12). In the absence of this solution, a new procedure 
for rotating the design parameter space must be defined. The procedure is 
first explained for the case of a two-design-paranieter problem, and then it is 
extended to the general multi-design-parameter problem. 

Figure 1 shows the design parameter space for a two-desi gn-parameter 
problem. In the figure* the constraint function values f^, f 7 ^ are defined 

as follows: 


f o = §"> 


t 


t 




( 


i 


r- il r- / . y II U\ 

fj^ _ ^(f E i 1 > Sj) 




£ 2 * ^ ~ ^ ~ 2 9 ^ 2 ^ 

In the modified procedure, the chord method, used in Equation (7) to satisfy the 
constraint condition, is used to rotate the design parameter space. The rota- 
tion angle 66™ +1 given by 


66 


n+1 

M 


reef* - f^) 

-1 l o 

tan 

e 


(14) 


is used to rotate the coordinate system, where the subscript M indicates that 
the modified scheme is used. The angle 60^ is now compared to the correspond- 
ing rotation angle 60 n+ l used in the original scheme and given by 


60 


n+1 


tan 


-1 



€ 



(15) 
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This comparison shows that the tern f™ - f™ in the original schene is replaced 
by e/C in the modified schene. Therefore, the modified schene may be viewed as 
the original schene with the exception that the exact value for i^ is replaced 
by an approximate estimate in which the gradient of f in the direction of the 
axis, is not calculated but is estimated using the same proportionality 
constant used in the chord method of Equation (7). Thus, 



This is applicable for both the two-design-parameter problem and the general 
multi-design-paraneter problem. 

In the optimization schene developed here, corrective increments are 
applied to the design parameter solutions every few iterations of updating the 
flow solutions. For convergence to occur, the signs of the increments must be 
chosen correctly to allow the iterative solution to approach the desired 
solution. The magnitudes of the increments are dependent on the computational 
constants c^, C 2 , and C. Because the design parameters are updated 
frequently during the iterative process, we are not concerned with determining 
the incremental step sizes that lead to the highest short-term convergence 
rate. In fact, this may be difficult to define, since the flow variable 
solutions are continuously changing during the iterative process. Our aim is 
to achieve design parameter convergence over a long tern defined by the number 
of iterations required for the flow solution convergence. A wide range of 
incremental step sizes should produce the desired convergence properties over 
many iterations, even though convergence properties over a few iterations may 
differ. These comments apply to both of the schemes described above for deter- 
mining the design parameter space rotation. The direct procedure for determin- 
ing the design parameter space rotation in the original scheme is replaced by 
an iterative procedure in the modified scheme. Since this rotation is updated 
frequently during the iterative process, this replacement should have no sub- 
stantial effect on the overall convergence of the solution. 

A potential problem exists when the modified scheme is used for rotating 
the design parameter axes. This problem is now discussed and suggestions for 
overcoming it are then presented. 

In the first aN- 1 iterative steps of solving the problem, the coordinate 
system in the design parameter space coincides with the original unrotated 
design parameter space P^, P 2 , . . ., P^. At the aN^ iterative step, a new 


4 


I 


i 


I 


4 


4 


4 


4 
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rotated coordinate system is determined. When Equation (13) for determining 

is used, we are guaranteed that the vector i^ points in the direction in 

which the constraint function increases. Consequently, the use of Equation (7) 

will cause the iterative solution to approach the constraint surface. When 

Equation (13) is replaced by Equation (16) for determining G^ \ there is a 

possibility that the computed vector i^ points in the direction in which the 

constraint function decreases. In this case, the assumption that C is positive 

is wrong, and using it will cause the solution to diverge. This occurs if the 

aN-1 

vector e^ is nearly in the direction of -Vf * . That is, if the quantity 

' Ivf 4 "- 1 ! ' §1 

is close to unity. The probability of this occurring is approximately 1:4 in a 
two-design-paraneter problem and is reduced further as the number of design 
parameters increases. There are two suggested approaches for overcoming this 
problem. In the first approach, the initial few iterations are performed 
using the original scheme for determining G™ by Equation (13) in order to 
determine the correct initial directions for the axis. This may then be up- 
dated using the modified scheme, Equation (16), in the rest of the computation. 
Realizing that the probability for the potential problem to occur is small, 
the second approach uses the modified scheme from the beginning of the 
computation. If divergence does occur, then the constraint function is 
redefined to be equal to the negative of the original constraint function and 
the problem is solved again. 
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OPTIMIZATION PROGRAM 


The optimization scheme described above is applicable to general 
aerodynamic problems, and can be used in conjunction with different analysis 
codes. The scheme was tested by applying it to the problem of optimizing 
advanced propeller designs.^ In these tests it was used in conjunction with 
code NASPROP-E which is an analysis code that computes the flow around an 
advanced propeller. A description of the scheme used in solving the flow 
equations around the propeller is given in Reference 6. 

In the following sections we will focus on describing the elements of the 
optimization program. It will be necessary to make some reference to code 
NASPROP-E; however, this reference will be minimized. Names of subroutines or 
COMMON blocks belonging to NASPROP-E will be bracketed to indicate that they 
do not belong to the optimization program and that they are replaceable by 
other subroutines or COMMON blocks when the optimization scheme is used in 
conjunction with other analysis codes. A manual describing the main elements 
of code NASPROP-E is given in Reference 8. 

Flow Chart 

In Figure 2, a brief flow chart describing the overall flow of the program 
is presented. 

Program Input 

The input required by the optimization program is described in this 
section. A dictionary of the input variables is provided, followed by a 
description of the input data format and an example of a set of input data. 

Dictionary of Input Variables 

The variables which are input to the optimization program are described in 
the following list. 

BARMl(L) , L-l ,NPARM 

Minimum allowable value for the Lth design parameter, in rotated 
space, during the iterative process. 

BARM2(L) ,L=1,NPARM 

Maximum allowable value for the Lth design parameter, in rotated 
space, during the iterative process. 
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Figure 2. Program Flow Chart 
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Cl 

C2 

CAA 

CNSTR 

DELPRM 

DPARM(L) 

DPARMX(L) 


DSIN(M) 


ISLCT(L) 

IRSC 

=0 

=1 

ISVC 

= 0 
=1 

H STEP 
NDSIN 


NHBL12 =1 
=2 

NPARM 

=0 

=1 

>1 


Incrementing factor for the optimization scheme. 

Decrementing factor for the optimization scheme. 

Proportionality constant for the chord method. 

Value of constraint function. 

Small positive Incremental value used to perturb the design 
parameters. 

,L=2, NPARM 

Initial incremental step size for the Lth design parameter. 

,L=1, NPARM 

Maximum allowable incremental step size for the Lth design 
parameter in rotated space. 

,M**1 , NDSIN 

Initial values for the geometrical parameters. These geometrical 
parameters are coefficients of the shape functions defined in 
subroutine OPTSHP. 

,L=1, NPARM 

The Lth component of the design parameter vector is equal to the 
ISLCT(L) component of the vector of geometrical parameters. 

This is a parameter which identifies the type of restart data. 

The restart data was saved from a regular analysis run. 

The restart data was saved from an optimization run. 

This is a parameter that identifies the type of data being saved. 
The present run is a regular analysis run. 

The present run is an optimization run. 

Number of iterative steps at which the design parameters are 
periodically updated. 

Number of geometrical parameters used to perturb the aerodynamic 
shape from the baseline configuration. A number, NPARM, of the 
geometrical parameters are allowed to vary. These are the design 
parameters. The rest of the geometrical parameters are fixed. 
NDSIN must have a value of 1 or greater. 

The hub shape is optimized 
The blade shape is optimized 
Number of design parameters 
Regular analysis problem 
Constraint problem 
Constrained optimization problem 
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NPRNT Data associated with the iterative history of the solution is 

printed every NPRNT iterative steps. 

NTRNOP =0 Coordinates will not be rotated. NPARM problems are solved in 
parallel. This option should be used only if the constraint 
surface is a planar surface and the first design parameter 
coordinate is normal to that surface and increases in the 
direction of increasing constraint function. In general, this 
option is not useful. 

=**1 Coordinates will be rotated. NPARM problems are solved in 
parallel. is determined by Equation (16). 

“1 Coordiates will be rotated. NPARM + 1 problems are solved in 
parallel. is determined by Equation (13). 

RRRANG Parameter used in the definition of the shape functions given in 

subroutine OPTSHP. 

RRRC Parameter used in the definition of the shape functions given in 

subroutine OPTSHP. 


Input Data Format 

The input data for the optimization program is read from subroutine 
OPTRED. The following reading sequence is used with the format given below. 
READ (55,100) 

READ (55,200) NPARM, MSTEP, NTRNOP, NDSIN, NPRNT, IRSC, ISVC 
READ (55,100) 

READ (55,200) (ISLCT (I), 1=1, MMPARM) 

READ (55,100) 

READ (55,300) (DSIN (I), 1=1, NDSIN) 

READ (55,100) 

DO 1 N=1 , MMPARM 

READ (55,300) DPARM(N) , DPARMX(N) , BARMl(N), BARM2(N) 

1 CONTINUE 
READ (55,100) 

READ (55,300) CNSTR 
READ (55,100) 

READ (55,300) CAA, Cl,C2, DELPRM 
READ (55,100) 

READ (55,300) RRRC, RRRANG 
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READ (55,100) 

READ (55,200) NHBL12 
100 FORMAT (IX) 

200 FORMAT (8110) 

300 FORMAT (8F10.0) 

Note: MMPARM = 1 If NPARM = 0 and MMPARM = NPARM otherwise. When NPARM = 

0, the variables ISLCT(l), DPARM(l), DPARMX(l), 3ARM1(1), and BARM2 ( 1 ) must be 
specified in the input data, even though they are not used in this particular 
case. The variable DPARM(l) is always not used; however, a value must be 
specified for it in the input data. 

An additional set of data is read from the subroutines of the analysis 
code used in conjunction with the optimization program. This data includes 
the baseline description of the aerodynamic configuration, flow conditions 
(for example, the free stream Mach number value), computational parameters 
associated with the analysis code, and the desired number of iterations, 

NMAX. For the propeller problem, this set of data is decribed in Reference 3. 

Input Data Example 

The following is an example of a set of input data: 


NPARM 

MSTEP 

NTRNOP 

NDSIN 

NPRNT 

IRSC 

ISVC 

2 

40 

-1 

3 

20 

0 

0 

ISLCT(l) 

ISLCT( 2) 






1 

3 






DSIN(l) 

DSIN( 2) 

DS IN ( 3 ) 





2.0 

0.0 

2.0 





DPARM(I) 

DPARMX(I) 

BARMI(I) 

BARM2(I) 




0.0 

1.0 

-50.0 

50.0 




0.5 

1.0 

-20.0 

20.0 




CNSTR 







1.7 







CAA 

Cl 

C2 

DELPRM 




3.0 

1.2 

0.6 

0.0001 




RRRC 

RRRANG 






0.25 

0.5 






NHBL12 







2 
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COMMON Blocks 


The COMMON blocks used in the optimization program are given below and the 
variables appearing in these blocks are defined. This is followed by a COMMON 
block/subroutine cross reference. 

Dictionary of COMMON Block Variables 

COMMON BLOCK PARM1 

C0MM0N/PARM1/NP, NB, MM, Ml, M2, RADI 

The variables belonging to this COMMON block are variables used in code 
NASPROP-E. 

COMMON BLOCK PARM4 
COMMON / P ARM 4 /N PRN T 

NPRNT Output and plotting data is written out every NRPNT iterative 

steps . 

COMMON BLOCK PARM5 

C0MM0N/PARM5/DBET0(I, 2) , DBTINC(I,L), DSIN(L), NDSIN, ISLCT(L), RNO(II) , 
DRNINC(II,L) 

All dimensions, I, appearing in this COMMON block should have a value of 
NB or greater, where NB is the number of planes In the radial direction 
defining the blade geometry. All dimensions, II, appearing in this COMMON 
block should have a value of NN or greater, where NN is the number of points 
in the axial direction defining the nacelle geometry. All dimensions, L, 
appearing in the COMMON block should have a value of NDSIN or greater. 


DBET0( J, 2) 

Input or baseline blade twist at each radial station 
relative to the 75% span station. 

DBTINC( J,K) 

Perturbation blade twist at each radial station, 
corresponding to the Kth geometrical shape function. 

DRNINC( J,K) 

Perturbation nacelle radius at each axial station, 
corresponding to the Kth geometrical shape function. 

DSIN(J) 

Jth geometrical parameter. 

NDSIN 

Number of geometrical parameters. 

ISLCT(J) 

The Jth component of the design parameter vector is 
equal to the ISLCT(J) component of the vector of 
geometrical parameters. 

RNO(J) 

Input or baseline nacelle radius at each axial station. 
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COMMON BLOCK PARM6 
C0MM0N/PARM6/RRR(I), RRRC, RRRANG 

The dimension I in this COMMON block should have a value of NB or greater, ^ 

where NB is the number of planes in the radial direction defining the blade 
geometry. 

RRR(J) Radius of radial stations on the blade. 

RRRANG Parameter used in the definition of the shape functions. ^ 

RRRC Parameter used in the definition of the shape functions. 

COMMON BLOCK PARM7 

C0MM0N/PARM7/IRSC, ISVC, NNCNT * 

IRSC Parameter defining the type of restart data, 

ISVC Parameter defining the type of data being saved. 

NNCNT Counter 

COMMON BLOCK PARM40 ( 

COMMON / PARM40 /NHBL12 

NHBL12 =1 The hub shape is optimized 

=2 The blade shape is optimized 

< 

COMMON BLOCK PAR171 

C0MM0N/PAR171/PARM( I) , DPARM(I), EP(I), FP(l), BARMl(I), BARM2( I) , DPARMX(I), 

XNRM(I), TRAN (1,1), TRANT(I,I) , TRAN0(I,I), NPARM , MSTEP, 

CAA1, CAA2, CAA, DELPRM , EO, FO, PARMOC(I), NTRNOP 
All dimensions, I, appearing in this COMMON block should have a value of 
NPARM or greater. 

BARMl(L) , L=1 , NPARM 

Minimum allowable value for the Lth design parameter, in 
rotated space, during the iterative process. 

BARM2CL) ,L-1, NPARM 

Maximum allowable value for the Lth design parameter, in 
rotated space, during the iterative process. 

CAA Proportionality constant for the chord method. 

CAA1 Incrementing factor for the optimization scheme. 

CAA2 Decrementing factor for the optimization scheme. 

DELPRM Small positive incremental value used to perturb the design ' 

parameters. 
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DPARM(L) 


DPARMX(L) 


EO 

EP(L) 


FO 

FP(L) 


MSTEP 


NPARM 

NTRNOP 

PARM(L) 

PARMOC(L) 


THAN 


TRANT 


TRANO 
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,L-1, NPARM 

Incremental step size for the Lth design parameter in rotated 
space . 

, L=1 , NPARM 

Maximum allowable incremental step size for the Lth design 
parameter in rotated space. 

Objective function corresponding to main solution. 

,L=1, NPARM 

Objective function corresponding to the Lth perturbed design 
parameter in rotated space. 

Constraint function corresponding to the Lth perturbed design 
parameter in rotated space. 

,L=1 , NPARM 

Constraint function corresponding to the Lth perturbed design 
parameter in rotated space. 

Number of iterative steps at which the design parameters are 
periodically updated. 

Number of design parameters. 

Parameter defining scheme options 
,L=1, NPARM 

Lth design parameter in rotated space. 

,L=1, NPARM 

Lth design parameter in original unrotated space. 

The first NPARM x NPARM elements of this two-dimensional array 
form the transformation matrix which operates on a vector with 
components defined relative to the rotated coordinate system of 
the present iteration to give the corresponding components rela 
tive to the rotated coordinate system of the previous iteration 
The first NPARM x NPARM elements of this two-dimensional array 
form the transpose of the transformation matrix defined by the 
elements of TRAN. 

The first NPARM x NPARM elements of this two-dimensional array 
form the transformation matrix which operates on a vector with 
components defined relative to the present rotated coordinate 
system to give the corresponding components relative to the 
original unrotated coordinate systef. 
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COMMON BLOCK PAR172 
COMMON/PAR172/CNSTR 

CNSTR Constraint function. 

COMMON BLOCK PAR176 
COMMON /PAR17 6 /LLMIN, LIMAX 

LLMIN Number of initial iterative step. 

LLMAX Number of final iterative step. 

COMMON BLOCKS (BASE), (COUNT), (CRIDC) 

These COMMON blocks are among the NASPROP-E COMMON blocks. 

COMMON BLOCK PWEF 
COMMON/PWEF/CPI, ETA, II 
CPI Power 

ETA Efficiency 

II Counter giving the number of times output data has been 

written . 

COMMON BLOCK TRAC 

COMMON /TRAC /TIME( I), RESMX(I), RESAV(I), POW(I), EFF(I) , DSNPR(I,J) 

All dimensions, I, appearing in this COMMON block should have a value of 
KK or greater, where 

KK = NMAX / NPRNT , 

NMAX is the desired number of iterations, and the output data is written out 
once every NPRNT iterations. The dimension, J, appearing in the COMMON block 
should have a value of MMPARM or greater where MMPARM = max(l ,NPARM) . 


TIME(K) 

RESMX(K) 

RESAV(K) 

POW(K) 


,K=1 , KK 

Number of iterative steps for the Kth data point. 
,K=1 , KK 

Maximum residual at Kth data point. 

,K=1 , KK 

Residual Euclidean norm at Kth data point. 

,K=1, KK 

Power coefficient at the Kth data point. 
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EFF(K) 


»K=1 , KK 

Efficiency at the Kth data point. 

DSMPR(K, L) ,K=1 , KK ,L=1, NPARM 

Iterative solution for the Lth design parameter at the Kth data 
point . 


COMMON BLOCK VAR1 

This COMMON block is used to save memory storage requirements. Variables 
contained in this COMMON block in subroutine (MAIN4) are associated with mesh 
generation. These variables are replaced with flow field solutions in all 
other subroutines which contain the COMMON block VAR1. In subroutine (MAIN4) 
the COMMON block appears as 

C0MM0N/VAR1/GRIDX, GRIDY, GRIDZ, GRIDR, GRIDP 
where GRIDX, GRIDY, GRIDZ, GRIDR, and GRIDP are arrays that contain mesh infor- 
mation. In all other subroutines containing COMMON block VARl, it appears as 
COMMON / VARl / Q , S, QM1 , SCRAT 

where Q, S, and QM1 are arrays that contain flow field solutions, while the 
array SCRAT is not used. 


COMMON Block/Subroutine Cross Reference 

main CORE COREIH ERVAL FRVAL OPTGOM OPTITR OPTRED OPTSHP PARMS SWTCH1 SWTCH2 


PARM 1 X 

PARM4 
PARM 5 X 

PARM6 X 

PARM7 X 

PARM 40 

PARI 71 X X X X 
PAR172 X 

PARI 7 6 X 

( BASE) X 

(COUNT) 

( CRIDC) 

PWEF X XX 

TRAC X 

VARl X X 


X 


X X 

X 

X X 

X X 

X 
X 

X 

X 


X X 
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Program Subroutines 
Subroutine Cross Reference 

For each subroutine, the following table lists the subroutines that it 
calls and the subroutines that call it. 


Subroutine 
Main Program 


CORE 

C0RE1H 

ERVAL 

FRVAL 

GRAM 

OPTGOM 

OPTITR 


OPTRED 

OPTSHP 

PARMS 

SWTCH1 

SWTCH2 

TRNSFR 

XMORM 

(BNDRY2) 

(FORCE) 

(GRID) 

( INPUT) 
(METRIC) 
(PRMESH) 
( STEP) 


Called By 

Calls 

CORE, C0RE1H, OPTITR, 
OPTRED, PARMS, SWTCH1, 
SWTCH2 

Main Program 

GRAM, XMORM 

Main Program 

TRNSFR 

OPTITR 

(FORCE) 

OPTITR 

CORE 

(FORCE) 

OPTITR 

(GRID), (INPUT), 
(METRIC), (PRMESH) 

Main Program 

ERVAL, FRVAL, OPTGOM, 
PARMS, SWTCH1, SWTCH2, 
TRNSFR, (BNDRY2), (STEP 


Main Program 
MREAD 

Main Program, OPTITR 
Main Program, OPTITR 
Main Program, OPTITR 
C0RE1H, OPTITR 
CORE 
OPTITR 

ERVAL, FRVAL 

OPTGOM 

OPTGOM 

OPTGOM 

OPTGOM 

OPTITR 
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Subroutine Descriptions 

In this section, a brief description of each of the subroutines of the 
optimization program Is given. Also described briefly are the subroutines of 
Program NASPROP-E that were referenced In this report. 

SUBROUTINE CORE 

This subroutine computes the transformation matrix TRAN and Its transpose 
TRANT. These two matrices relate the components of a vector relative to the 
current rotated coordinate system to Its components relative to the rotated 
coordinate system of the previous iteration. 

SUBROUTINE COREIH 

This subroutine computes the new iterative solution for the design 
parameters. All the dimensions of the variables appearing in the DIMENSION 
statement in this subroutine should be given a value equal to or greater than 
NPARM, the number of design parameters. 

SUBROUTINE ERVAL 

This subroutine computes the value of the objective function. The user is 
required to add neccessary COMMON blocks and necessary calls to other 
subroutines which allow the computation of the objective function. 

Argument List 

E (Output) Objective function 

SUBROUTINE FRVAL 

This subroutine computes the value of the constraint function. The user 
is required to add necessary COMMON blocks and necessary calls to other 
subroutines which allow the computation of the constraint function. 

Argument List 

F (Output) Constraint function 

SUBROUTINE GRAM 

This subroutine uses the Gram-Schmidt orthogonalization process to con- 
struct a set of orthonorraal vectors from a set of linearly independent vectors. 
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Argument List 

E (Input) Two-dimensional array. Set the first dinension equal to 

the first dinension of TRAMT in COMMON block PAR171. Set the 
second dinension equal to 1. The first N columns of array E 
contain in their N first rows the comnonents of the innut set of 
linearly independent vectors. 

ET (Outnut) Two-dimensional array. Set the first dinension equal to 

the first dinension of TRAN in COMMON block PANT 7 1 . Set the 
second dinension equal to 1. The first N columns of array ET 
contain in their N first rows the components of the output set of 
orthonormal vectors. 

N (Input) Number of vectors. 

NBIG (Input) The comnonent number with the biggest magnitude in the 
first column of the input array E. 

SUBROUTINE OPTGOM 

This subroutine Is called so that the comnutati onal mesh may be 
generated. Subroutine OUTGO' \ c n 1. 1 ?> subroutine (PT!ES.:0 , •./hi ch does the act ‘ 
computation of t Vis? mesh. The storage area allocated for COMMON block VARl is 
used for mesh computations. Array QU1 , appearing as a variable in COMMON 
block VARl , stores the flow solution. For that reason, array QM1 is written 
to tape 60 before subroutine (PRMESH) is called. Following the call to 
subroutine (PRMESH), the flow solution is restored in array QM1 by reading it 
from tape 60. 

Argument List 

LLKK (Input) 

=1 the mesh is computed 

=2 mesh dependent variables are computed 
SUBROUTINE OPTITR 

This subroutine updates the main flow solution or the perturbed flow 
solution by performing MSTEP iterative steps for fixed values of the design 
parameters. The dimensions of the variables appearing in the DIMENSION 
statement in the subroutine should be given a value equal to or greater than 
NPARM, the number of design parameters. 






( 


{ 
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Argument List 

N ( Input ) 

=0 Main solution is updated 

i0 Perturbed solution is updated 

NCYCB (Input) Number of cycles for updating the design parameters. 

SUBROUTINE OPTRED 

This subroutine reads the input data associated with the optimization 
scheme. The data is read from tape 55. 

SUBROUTINE OPTSHP 

This subroutine defines the shape function used to perturb the baseline 
aerodynamic configuration. The shape functions are defined by the 
two-dimensional array DBT1NC in the case of optimizing the blade angle 
distribution and by the two dimensional array DRNINC in the case of optimizin 
the hub shape. The first argument of these arrays defines a radial position 
along the blade while the second argument defines a particular shape function 


SUBROUTINE PARMS 

This subroutine sets the design parameters, which are contained in the 
elements of the array PARK, equal to the geometrical parameters, contained in 
the elements of the array DSIN, if NPL1R2 =1. If NPL1R2 = 2, the subroutine 
sets the geometrical parameters equal to the design parameters. 


Argument List 


NPL1R2 


PARM 

NPARM 


(Input ) 

=1 Elements of PARM are set equal to elements of DSIN 
=2 Elements of DSIN are set equal to elements of PARM 
(Input or Output) 

Vector of design parameters 
(Input ) 

Number of design parameters 
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SUBROUTINE SWT CHI 

The optimization scheme requires that a number of solutions for different 
flow problems be obtained in parallel. The argument N refers to the main 
solution if N=0, and to the Nth perturbed solution if N^O. The iterative 
solutions for these different problems are stored in the array QM1. Prior to 
updating a particular solution, it is temporarily stored in array Q. The code 
then updates array Q, replacing the old solution by the solution at the new 
iterative step. This solution is again stored in array QM1. Subroutine 
SWTCH1 stores the solution corresponding to the Nth problem in array Q prior 
to updating it. 

SUBROUTINE SWTCH2 

This subroutine stores the solution corresponding to the Nth problem in 
array QMl after it is updated. 

SUBROUTINE TRNSFR 

A transformation matrix relating two coordinate systems, (1) and (2), is 
used in this subroutine to compute the vector components relative to 
coordinate system (2) from the vector components relative to coordinate system 
( 1 ). 

Argument List 

VIN (Input) Vector whose components are defined relative to coordinate 
system 1. 

VOUT (Output) Vector whose components are defined relative to 
coordinate system 2. 

T (Input) Transformation matrix. The first dimension in this array 

should be given a value equal to the first dimension of the 
variable TRAN in COMMON block PAR171. 

N (Input) Dimension of the space under consideration. 

SUBROUTINE XNORM 

This subroutine computes the components of the unit vector normal to a 
constraint surface contour. 


TR-444/07-88 


24 



Argument List 

DP (Input) Small positive number used as a perturbing parameter. 

FPNP (Input) Vector whose elements contain the values of the constraint 

function at positions perturbed from the point of interest by DP 
along each of the coordinate directions. 

FONP (Input) Value of the constraint function at the point of interest. 
X (Output) Unit vector normal to the constraint surface. 

N (Input) Number of design parameters. 

CAA (Input) Proportionality constant used in the chord method. 

NTRNOP (Input) Parameter determining the procedure for computing the unit 
normal vector. 

SUBROUTINE (BN DRY 2) 

This subroutine satisfies the flow boundary conditions. 

SUBROUTINE (FORCE) 

This subroutine computes the power and efficiency. 

SUBROUTINE (GRID) 

This subroutine reads in the computational mesh from TAPE11. 

SUBROUTINE (INPUT) 

This subroutine reads the input data for the aerodynamic analysis code. 
SUBROUTINE (METRIC) 

Quantities which are dependent on the computational mesh are computed in 
this subroutine. 

SUBROUTINE ( PRMESH) 

This is the main subroutine for the mesh generation program. Its function 
is to call other subroutines which compute the mesh. The mesh data is 
written on TAPE11. 

SUBROUTINE (STEP) 

This is the subroutine which computes the flow solution at a new iterative 
step. 
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Program Output 


An iterative history showing the iterative time step, the residual, the 
power, the efficiency, and the design parameters is printed out. Similar dat 
is saved on TAPE40 for the purpose of producing line plots. This data is 
written in Subroutine (OUTPUT). 
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APPENDIX 
LIST OF SYMBOLS 


= incrementing factor for optimization scheme [see Equation (9)] 

= decrementing factor for optimization scheme [see Equation (9)] 

= positive constant for chord method [see Equation (7)] 

= power coefficient 
= desired power coefficient 

= unit vector along the P axis 

£ 

= objective function 
= constraint function 

= solution of the flow governing equations 

= component of Vf relative to rotated coordinate system 
= unit vector along the P^ axis with components defined relative to the 
unrotated design parameter coordinate system 
= unit vector along the P^ axis with components defined relative to the 
rotated design parameter coordinate system 
= number of design parameters 
= vector of design parameters 

= vector of design parameters relative to rotated coordinate system 
t Vi 

= a, component of design parameter vector 
t-ll 

= £ component of design parameter vector relative to rotated coordinate 


system 
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= incremental vector used to update the vector of design parameters 
<5P 3 maximum incremental value allowed in updating the design parameters 

H3 X 

AN = number of iterative steps at which P is periodically updated 

e “ small positive incremental value used to perturb the design 

parameters 
n = efficiency 

■ flow iterative solution 
Superscripts 

n “ iteration number 

* =* optimum value 

= rotated coordinate system 

Subscripts 

M “ coordinate system rotated by the modified scheme 
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